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ABSTRACT 


An  axisymmetric  boundary  element  flow  algorithm  is  coupled  with  a  finite 
element  structural  analyzer  to  perform  interactive  calculations  of  the  growth 
and  subsequent  collapse  of  an  explosion  bubble  near  a  submerged  compliant 
structure.  The  validity  of  the  program  is  established  through  direct  compar¬ 
isons  with  published  experimental  data.  A  parametric  study  of  the  interaction 
between  a  growing  and  collapsing  bubble  and  a  spherical  shell  is  presented. 
The  results  show  that  if  the  stiffness  of  the  shell  is  sufficiently  low,  the  mass 
of  the  shell  is  a  critical  parameter  in  the  collapse  problem.  If  the  mass  of 
the  structure  is  high,  a  reentrant  jet  forms  and  is  directed  towards  the  shell. 
As  the  mass  of  the  spherical  shell  is  decreased,  the  collapse  becomes  spherical 
with  no  jet  formation.  At  the  lowest  structural  mass  for  which  calculations 
are  performed,  a  jet  directed  away  from  the  structure  begins  to  form.  The 
ratio  of  the  depth  of  submergence  to  bubble  maximum  radius  was  also  found 
to  be  a  critical  parameter  in  the  collapse  problem.  When  this  ratio  is  large 
(greater  than  100),  the  collapse  is  driven  by  interaction  forces.  However,  for 
shallow  submergence,  buoyancy  effects  become  more  important  than  interac¬ 
tion  forces. 

ADMINISTRATIVE  INFORMATION 

This  work  was  funded  by  the  David  Taylor  Research  Center  Independent 
Research  (IR)  program  under  Work  Unit  1750-161  for  fiscal  year  1991, 

INTRODUCTION 


BACKGROUND 


The  interaction  between  a  growing  and  collapsing  bubble  and  a  sub¬ 
merged  compliant  structure  is  relevant  to  a  number  of  engineering  prob¬ 
lems.  Important  examples  include  cavitation  bubbles  in  fluid  machinery  and 
pulsating,  collapsing  underwater  explosion  bubbles  near  naval  vessels.  Ex¬ 
perimental  and  theoretical  studies  have  shown  cavitation  bubble  collapse 
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to  be  characterized  by  extremely  high  local  velocity  and  overpressure,  and 
high-speed  liquid  jet  formation  (Roberson  and  Crowe1).  Recent  studies  have 
shown  bubble  centroid  motion  as  well  as  jet  strength  and  direction  to  be 
greatly  influenced  by  the  presence  and  nature  of  any  nearby  boundaries.  The 
jet  is  known  to  cause  surface  erosion  and  pitting  in  turbomachinery.  Ship 
propellers  and  other  propulsion  devices  are  especially  susceptible  to  this  type 
of  damage  (Arndt2).  Collapsing  cavitation  bubbles  are  also  responsible  for 
structural  vibration,  noise  ^nd  degradation  of  operating  efficiency. 

The  detonation  of  a  submerged  explosive  is  characterized  by  the  conver¬ 
sion  of  the  original  explosive  material  into  a  gas  sphere  at  extremely  high 
temperature  and  pressure.  The  inertia  of  the  water  surrounding  the  bubble 
and  the  compressible  explosion  products  form  a  damped  oscillatory  system 
(Cole3).  This  oscillating  bubble  and  entrained  fluid  can  induce  low  frequency 
beam- like  flexural  motions  in  a  nearby  ship  or  submarine  (Hicks4,5).  This 
whipping  phenomenon  can  lethally  damage  both  ship  structure  and  onboard 
equipment.  As  was  the  case  for  the  cavitation  bubble,  the  explosion  bubble 
is  greatly  affected  by  the  presence  of  any  nearby  boundaries.  Under  certain 
conditions,  explosion  bubble  collapse  is  also  accompanied  by  the  formation  of 
a  high  speed  reentrant  jet  (Snay6).  Here  again,  the  jet  strength  and  direction 
depend  upon  the  proximity  and  compliance  of  any  nearby  boundaries.  Stud¬ 
ies  on  the  damaging  potential  of  underwater  explosions  have  often  focused 
on  the  initial  shock  wave  emitted  by  the  explosion.  However,  many  investi¬ 
gations  have  shown  that  bubble  collapse  damage  can  be  more  extensive  than 
shockwave  damage. 

PRIOR  WORK 

To  study  the  fundamental  aspects  of  these  problems,  a  number  of  theo¬ 
retical,  numerical,  and  experimental  investigations  have  been  conducted  to 
examined  the  growth  and  collapse  of  vapor  bubbles  near  plane  boundaries. 
Collapse  near  a  rigid  wall  is  often  characterized  by  the  formation  of  a  reen¬ 
trant  jet  directed  toward  the  wall  (Benjamin  and  Ellis7).  When  the  plane 
boundary  is  a  free  surface,  a  reentrant  jet  is  also  formed,  but  in  this  case, 
the  jet  is  directed  away  from  the  boundary  (Blake  and  Gibson8). 
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This  complex  interaction  between  the  collapsing  bubble  and  nearby  bound¬ 
aries  has  led  a  number  of  investigators  to  suggest  application  of  a  relatively 
flexible  surface  coating  in  order  to  redirect  the  reentrant  jet  and  hence  re¬ 
duce  bubble  damage.  Shima,  Tomita,  Gibson,  and  Blake9  photographed 
the  growth  and  collapse  of  spark  generated  vapor  bubbles  near  planar  sur¬ 
faces  covered  with  elastomeric  coatings.  Their  study  indicated  that  for  some 
coating  properties  the  jet  was  indeed  redirected  away  from  the  solid  bound¬ 
ary.  Through  this  same  study,  detailed  histories  of  the  motion  and  shape  of 
the  bubbles  are  provided.  Numerical  simulations  of  the  interaction  between 
a  growing  and  collapsing  vapor  bubble  and  a  planar  compliant  surface  have 
been  conducted  by  Duncan  and  Zhang10,11.  In  these  simulations,  a  boundary- 
integral  flow  algorithm  was  coupled  to  a  finite  difference  representation  of  a 
spring-backed  membrane.  By  equating  pressure  and  velocity  conditions  in 
the  flow  and  on  the  structure  at  the  fluid-structure  interface,  Duncan  and 
Zhang  achieved  a  completely  interactive  simulation.  The  results  of  these  sim¬ 
ulations  were  in  qualitative  agreement  with  the  experimental  work  of  Shima, 
Tomita,  Gibson  and  Blake9. 

In  the  present  study  Duncan  and  Zhang’s  numerical  method  is  extended 
to  include  an  axisymmetric  finite  element  structural  idealization.  The  finite 
element  model  uses  isoparametric  continuum  elements  and  a  linear-elastic 
material  model.  This  structural  model  allows  one  to  include  more  detailed 
information  about  the  structure  and  more  complex  structural  geometries 
than  the  idealized,  planar,  spring-backed  membrane  used  by  Duncan  and 
Zhang.  As  in  Duncan  and  Zhang’s  model,  inviscid,  incompressible,  potential 
flow  is  assumed  and  the  boundary  integral  method  is  again  employed  to  ef¬ 
fect  a  solution  of  the  potential  flow  problem.  By  matching  the  pressure  and 
motion  conditions  at  the  fluid-structure  interface,  an  interactive  simulation 
is  again  achieved.  In  addition,  the  present  work  includes  a  polytropic  model 
of  the  gas  inside  the  bubble  and  a  gravitational  field,  significant  parameters 
that  were  not  present  in  the  work  of  Duncan  and  Zhang.  These  extensions 
make  possible  simulations  of  the  growth  and  collapse  of  underwater  explosion 
bubbles  near  target  structures. 


3 


ORGANIZATION  OF  STUDY 

The  following  section  contains  an  account  of  the  theory  behind  the  bound¬ 
ary  integral  formulation  for  the  potential  flow  problem,  a  description  of  the 
finite  element  algorithm  used  to  compute  the  structural  response,  and  de¬ 
tails  of  the  numerical  implementation  for  the  bubble  growth  and  collapse 
problem.  Verification  calculations  are  then  given  of  the  collapse  of  a  bubble 
near  a  finite  element  idealization  of  a  spring  backed  membrane.  The  results 
of  these  calculations  are  compared  with  the  results  of  Duncan  and  Zhang. 
Next,  an  attempt  is  made  to  numerically  simulate  the  experiments  of  Shima, 
Tomita,  Gibson  and  Blake9.  A  parametric  study  of  the  growth  and  collapse 
of  an  explosion  bubble  near  a  compliant  sphere  is  then  presented.  Finally, 
the  results  of  the  present  work  and  suggested  directions  for  future  research 
are  summarized. 


THEORETICAL  MODEL 
MATHEMATICAL  FORMULATION 

A  schematic  showing  the  initial  configuration  of  the  bubble  and  a  spher¬ 
ical  structure  is  given  in  Fig.  1.  A  cylindrical  coordinate  system  is  used 
in  the  solution  of  the  interaction  problem.  The  radial  coordinate  is  r  and 
the  circumferential  angle  is  6.  The  problem  is  axisymmetric  about  the  z- 
axis.  The  spherical  structure  has  outer  radius  R,  and  is  centered  at  r  =  0, 
z  =  —  R„.  The  bubble  is  initially  (t  =  0)  spherical  with  radius  Rq,  and  its 
center  is  located  at  z  =  Zq.  The  pressure  in  the  fluid  far  from  the  bubble, 
Poo,  is  maintained  constant  during  the  calculation.  The  theoretical  model 
for  the  fluid  is  the  one  used  by  Duncan  and  Zhang.  For  clarity,  parts  of  the 
presentations  of  Duncan  and  Zhang  are  included  here  in  modified  form. 

The  fluid  motion  is  assumed  to  be  incompressible  and  inviscid  and  there¬ 
fore  satisfies  Laplace’s  equation: 


VV  =  0 


(1) 
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where  V2  is  the  Laplacian  operator  and  4>  is  the  Eulerian  velocity  potential. 
The  fluid  velocity  is  equal  to  the  gradient  of  the  velocity  potential,  u  =  V^. 
On  the  surface  of  the  bubble,  the  pressure  in  the  fluid  is  equal  to  the  pressure 
in  the  bubble,  Pg.  The  condition  imposed  on  4>  by  this  dynamic  uoundarv 
condition  can  be  written  as  Bernoulli’s  equation  in  material  derivative  form: 

^  =  jlW|!  +  +  «(*  -  Z.)  (2) 

where  Dj Dt  is  the  derivative  with  respect  to  time  following  a  fluid  particle, 
p  is  the  fluid  density,  and 
g  is  the  acceleration  due  to  gravity. 


The  noncondensible  gases  inside  the  bubble  are  assumed  to  have  a  polytropic 
behavior, 

PgV'1  =  const  (3) 

where  7  is  the  polytropic  constant  and  V  is  the  volume  of  the  bubble.  The 
kinematic  boundary  condition  on  the  surface  of  the  bubble  states  that  ma¬ 
terial  points  remain  on  the  surface  of  the  bubble: 

Dxc 

dT  =  (4) 

where  xc  is  the  position  vector  to  these  material  points. 

For  t  <  0,  the  structure  is  motionless  and  statically  loaded  by  uniform 
pressure,  Poo,  applied  by  the  fluid.  The  response  of  the  structure  to  the  fluid 
loading  is  computed  using  a  linear  elastic  finite  element  model.  The  struc¬ 
ture  is  represented  by  an  assemblage  of  discrete  elements,  interconnected  at 
nodal  points  on  the  element  boundaries.  Each  quadrilateral  element  has  four 
nodes  defined  in  orthogonal  two  dimensional  coord’  ate  space.  There  are  two 
displacement  degrees  of  freedom  (dofs)  associated  with  each  nodal  point.  An 
axisymmetric  element  formulation  is  used  such  that  each  element  represents 
the  portion  of  the  structure  that  is  swept  out  by  rotating  the  planar  element 
one  radian  about  the  z-axis  (see  Fig.  2).  The  solution  of  a  system  of  dynamic 
equations  of  motion  provides  the  displacement,  velocity,  and  acceleration  of 
each  nodal  point  in  the  element  assemblage  at  discrete  times.  The  system 
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equilibrium  equations  of  the  element  assemblage  are  derived  from  virtual 
work  considerations  and  are  of  the  form: 

Ml  +  A'£  =  R  (5) 

where  M  is  the  global  consistent  mass  matrix, 

K  is  the  global  stiffness  matrix, 

R  is  consistent  nodal  load  vector  obtained  from  the  total  pressure 
in  the  fluid  acting  on  the  structure, 

£  is  the  nodal  acceleration  vector  (two  components  per  node), 

£  is  the  nodal  displacement  vector. 


The  consistent  mass  matrix  contains  structural  mass  information  and  pro¬ 
vides  the  relationship  between  force  at  degree  of  freedom  i  and  acceleration 
at  dof  j.  The  global  stiffness  matrix  relates  the  applied  force  at  dof  i  to 
resulting  displacement  at  dof  j.  The  K  matrix  is  developed  from  material 
elastic  stress-strain  relations  and  linearized  strain-displacement  relations. 


During  the  collapse,  the  fluid  and  the  structure 


earized  equations  for  the  pressure  and  velocity  in  the  two  systems.  These 
equations  are  satisfied  at  the  undisturbed  position  of  the  fluid-structure  in¬ 
terface: 

=  J_±  (R\ 

dt  dn 

-Ps(x4,f)  =  -p~^  +  Poo  (7) 


where  (n  is  the  component  of  the  s  rface  displacement  in  the  direction 
normal  to  the  undisturbed  surface  of  the  structure, 
d£n/ dt  is  the  normal  velocity  of  the  structure  (positive 
outward  into  the  fluid),  and 

d(f>/dn  is  the  derivative  of  the  velocity  potential  in  the  direction 
normal  to  the  undisturbed  surface  of  the  structure. 


The  second  equation  is  the  linearized  Bernoulli  equation  and  Ps(x,,t)  is  the 
fluid  pressure  applied  to  the  structure  at  z,. 
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To  show  how  the  system  of  equations  can  be  advanced  in  time,  assume 
that  at  time  t,-  all  dependent  variables  are  known.  The  boundary  conditions 
on  the  surface  of  the  bubble,  Eqs.  2  and  4,  are  integrated  numerically  to  get 
the  position  of  the  surface  of  the  bubble  and  the  value  of  <f>  on  the  bubble  at 
time  £,+1.  The  finite  element  method  is  used  to  obtain  the  value  of  d(n/dt  = 
—d<t>/dn  on  the  surface  of  the  sphere  at  i1+i.  In  order  to  move  on  to  the 
next  time  step,  f,-+ 2,  the  values  of  V<j!>  must  be  known  on  the  bubble  surface 
for  use  in  Eqs.  2  and  4.  However,  at  this  point  only  the  value  of  d<f>/ds 
can  be  computed  (where  s  is  a  coordinate  along  the  bubble  surface).  Also, 
in  order  to  find  the  value  of  d<f>/dn  on  the  surface  of  the  structure  at  tl+2, 
the  pressure  must  be  known  on  the  surface  of  the  structure  at  t>+ 1  for  use 
in  the  finite  element  calculation  of  the  structural  response.  The  pressure  at 
U+ 1  can  be  obtained  from  Bernoulli’s  equation  (Eq.  7)  if  d<f>/dt  is  known. 
Thus,  the  value  of  <£  on  the  fluid-structure  interface  at  time  U+\  must  be 
found  to  obtain  a  finite  difference  approximation  for  dcfr/dt.  To  complete  the 
problem,  the  values  of  d<f>/dn  on  the  bubble  surface  and  <f>  on  the  surface  of 
the  structure  are  obtained  by  solving  Laplace’s  equation  in  the  form  of  an 
integral  equation  (Lamb12): 

/s.+s. 9{?’  ^ dS •  -  Ls.  <8> 

where  Sj,  is  the  surface  of  the  bubble, 

S,  is  the  interface  between  the  structure  and  the  fluid, 
p  is  a  field  point  that  is  on  the  surface  5  =  Si,  +  S„,, 
q  is  the  source  point  that  is  also  on  S,  g(p ,  q)  =  l/\p  —  <f|, 
n  is  the  normal  to  S  directed  outward  from  the  fluid,  and 
dSq  is  the  area  element  of  S  varying  the  point  q. 

Since  the  problem  is  axisymmetric,  the  positions  of  the  field  points  need 
only  be  considered  in  a  single  plane,  0  =  0.  Once  this  equation  is  solved,  the 
calculation  can  proceed  on  to  the  next  time  step  or,  with  a  companion  form 
of  this  integral  equation,  the  velocity  and  pressure  can  be  found  at  any  point 
in  the  fluid. 

A  special  treatment  is  needed  to  start  the  calculation  because  the  pressure 
distribution  on  the  structure  at  t  =  0  is  unknown.  The  following  procedure 
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was  used.  From  the  linearized  form  of  Bernoulli’s  equation  (Eq.  7)  evaluated 
on  the  structure  at  t  =  0,  it  can  be  seen  that  the  pressure  on  the  structure  is 
equal  to  the  partial  derivative  of  the  velocity  potential  with  respect  to  time 
evaluated  on  the  surface  of  the  structure.  The  calculation  of  the  initial  value 
of  d<f>/dt  on  the  structure  is  divided  into  two  steps.  In  the  first  step,  the 
initial  value  of  dfyfdn  on  the  bubble  surface  is  obtained  from  the  integral 
equation  (Eq.  8)  with  the  initial  value  of  4>  on  the  bubble  obtained  from  an 
integrated  form  of  the  Rayleigh- Plesset  equation  for  the  case  of  a  spherical 
bubble  growing  in  an  infinite  fluid  (see  Batchelor13), 


where  Pg o  is  the  pressure  inside  the  bubble  at  the  start  of  the  calculation  and 
a  =  Ro/Rmax ■  Since  the  structure  is  initially  motionless,  the  value  of  d<t>/dn 
on  the  structure  is  given  by 


d<t>  _  d(n 

dn  dt 


(10) 


The  values  of  dfyjdt  on  the  bubble  are  then  obtained  using  Bernoulli’s  equa¬ 
tion  in  the  form 

9t  2  \dn)  p  '  ’ 

where  the  fact  that  d<f>/ds  —  0  at  t  =  0  (s  is  the  direction  along  the  bubble 
surface)  has  been  used.  In  the  second  step  of  the  calculation,  it  is  noted 
that  the  partial  derivative  of  <f>  with  respect  to  time  also  satisfies  Laplace’s 
equation: 

vJ  (w)  “  °-  (12) 

Thus,  the  integral  equation  (Eq.  8)  is  valid  with  <f>  replaced  by  d<f>/dt.  On 
the  bubble,  d<i>jdt  is  obtained  from  the  first  step  described  above.  On  the 
surface  of  the  structure, 

d<l>  d  d<f> 

pm=ma^Tf  (I3) 

The  value  of  d<j>/dt  on  the  structure  and  thus,  the  pressure  on  the  structure 
at  t  =  0  can  be  obtained  by  solving  the  integral  equation  for  d<f>/dt. 
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NUMERICAL  IMPLEMENTATION 

In  the  numerical  model,  the  surface  of  the  bubble  is  approximated  by 
a  set  of  panels  each  of  which  is  obtained  by  rotating  a  straight  line  in  the 
0  —  0  plane  about  the  z-axis  (see  Fig.  3).  The  bubble  is  composed  of  nj  of 
these  panels.  The  interface  between  the  fluid  and  the  structure  is  modeled 
by  a  set  of  ns  equal  length  panels.  Field  points  (nodes)  are  taken  at  the 
positions  where  the  line  of  intersection  of  two  adjacent  panels  pierces  the 
0  =  0  plane.  A  predictor-corrector  scheme  sometimes  referred  to  as  Heun’s 
method  (see  Ferziger14)  is  used  to  perform  the  temporal  integration  of  the 
boundary  conditions  on  the  bubble.  The  r  and  z  coordinates  of  the  nodes  on 
the  bubble  and  the  corresponding  values  of  (p  at  t,+i  are  expressed  in  vector 
form  as: 

Predictor  step: 


1 

f  ri ) 

ui 

*«+i 

Pm  J 

>  -  i 

’  ■+■  (£«+i  —  U) ' 

w{ 

I|V#|2  +  ^ 

(14) 


Corrector  step: 


f  Wifi  1 

f  < 

'  =  {  zi 

l  Pm  J 

I 


P  + 

w\  +  wJi+1 

§|v^I2  +  |!v^+1|2  +  2^ 


(15) 


where  the  superscript  refers  to  the  nodes  and  the  subscript  refers  to  the  time 
step.  The  value  of  |V<£|  is  computed  from  the  derivative  of  <j>  in  the  direction 
normal  to  the  bubble  surface,  which,  in  turn,  is  obtained  from  the  solution 
of  the  integral  equation,  and  the  derivative  in  the  direction  tangent  to  the 
surface,  obtained  from  a  central  difference  scheme. 


In  the  numerical  solution  of  the  integral  equation  (Eq.  8),  the  values  of 
<f>  and  d<f>/dn  are  assumed  to  vary  linearly  in  the  0  =  0  plane  as  the  source 
point  is  varied  along  each  panel.  The  integral  equation  in  its  discrete  form 
can  be  written  as 
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-  £  (GNu’j<f>j  +  GN2'l-J<f>3+l) 


where  N  =  nc  +  nw  and 


.  .  rL3  JJ  _  p  /2jt  1 

Gu'3  =  /  —  /  -zr - z - < 

Jo  L 3  Jo  |p‘ —  gJ(r,  0)| 

rL>  .  13  f2ir  1 

G2’,,j  =  /  r(P)-—  /  -= - 4 - dOdP 

Jo  'D  Jo  |p._^(r,0)| 


-dddP 


-dOdP 


, . .  r13  Lj  -P  r2*  d  (  1  \ 

GNl™  =  /  r(P)———. —  /  —  — - = - )dddP 

Jo  L 3  Jo  dn  |p<  _  qj(r,  9)\J 

rU  13  r2ir  f)  /  1  \ 

GN 2'*’J  =  /  r(P)J-  -z-[— - = - )dddP  . 

Jo  v  ’ D  Jo  dn  \ |p*  —  qj(r,0)\/ 


In  these  equations,  the  length  of  panel  j  in  the  6  =  0  plane  is  given  by  L: ,  and 
P  is  the  distance  coordinate  along  the  panel.  The  parameter  a1  is  the  solid 
angle  within  the  fluid  subtended  by  the  fluid  surface  at  node  i.  The  integra¬ 
tions  in  the  0-direction  were  carried  out  analytically  following  the  method 
of  Jaswon  and  Symm15;  results  in  terms  of  elliptic  integrals  were  obtained. 
The  integrations  in  P  were  carried  out  numerically  using  Gauss-Legendre 
quadrature  techniques  for  the  regular  parts  of  the  integrals  and  the  quadra¬ 
ture  methods  of  Anderson16  for  the  singular  parts  of  the  integrals. 


The  structural  finite  element  program,  developed  by  Whang17,  has  been 
modified  to  include  isoparametric  continuum  elements  and  Newmark-Beta 
(see  Newmark  18)  time  integration  with  variable  time  step  size.  The  structure 
is  discretized  into  ne  quadrilateral  finite  elements,  interconnected  at  n,  nodal 
points.  Each  element  is  comprised  of  four  nodes.  For  each  structural  node  on 
the  fluid-structure  interface  there  is  a  corresponding  node  in  the  fluid  model. 
Coupling  of  the  fluid  and  structural  models  is  accomplished  by  exchanging 
pressure  and  velocity  information  between  these  pairs  of  fluid  and  structural 
nodes.  The  numerical  forms  of  the  boundary  conditions  at  the  flow-  structure 
interface  (Eqs.  6  and  7)  are 
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where  [^f]^  is  the  normal  velocity  of  structural  node  js,  and  is  the 

normal  derivative  of  the  velocity  potential  at  the  corresponding  fluid  node  jj. 
The  second  boundary  condition  merely  states  that  the  pressure  computed  at 
fluid  node  jj  is  applied  at  the  corresponding  structural  node  js. 


The  time  step  of  the  calculation  varied  during  each  run.  At  each  step, 
the  time  difference 

At  =  A<*>mai 

*  1  +  o.5<zLx 


was  computed,  where  A$max  is  a  constant  and  qmax  is  the  maximum  fluid 
velocity  on  the  surface  of  the  bubble  at  any  time  step.  A  minimum  time 
step,  A tmin,  and  a  maximum  time  step,  A 2mox,  where  chosen  such  that 
Atmin  =  AtmOx/200.  If  A tmin  <  Atv  <  A tmai,  then  A tv  was  used  as  the  time 
step.  If  A tv  >  A tmax  the  time  step  was  taken  as  Atmar,  while  if  A tv  <  A tm,n, 
the  time  step  was  taken  as  A tmj„. 


PROGRAM  VERIFICATION 

A  series  of  interactive  calculations  was  performed  to  assess  the  overall 
performance  of  the  coupled  flow  and  structural  algorithms.  This  was  ac¬ 
complished  by  making  direct  comparisons  between  the  results  of  the  present 
model,  employing  the  finite  element  method  and  the  results  of  Duncan  and 
Zhang’s  model,  using  the  finite  difference  method.  The  calculations  simu¬ 
lated  the  collapse  of  a  bubble  near  a  spring- backed  membrane. 

MODEL  DESCRIPTIONS 

A  schematic  showing  an  idealization  of  a  spring-backed  membrane  and  a  va¬ 
por  bubble  is  given  in  Fig.  4.  The  spherical  bubble  is  initially  at  rest  with 
radius  R0  and  is  centered  at  r  =  0,  z  =  Za.  The  pressure  inside  the  bubble, 
P0,  and  the  pressure  in  the  fluid  at  infinity,  P0 0,  are  held  constant  during  the 
calculation.  The  membrane  lies  in  the  r-9  plane  and  is  centered  at  r  =  0, 
z  —  0.  The  membrane  has  radius  Rm •  For  r  >  Rm  on  z  =  0,  the  boundary 
is  modeled  as  a  flat  rigid  wall.  A  constant  tensile  force,  T,  is  applied  to  the 
membrane  at  r  =  R^. 
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The  finite  difference  membrane  has  three  independent  variables,  mass  per 
unit  area,  m;  spring  constant  per  unit  area,  I\;  and  membrane  tension  per 
unit  length,  T.  The  independent  variables  for  the  flow  are  the  initial  radius  of 
the  bubble,  Ra,  the  initial  standoff  from  the  bubble  to  the  membrane,  Z0,  the 
pressure  diffe  rence,  AP  =  —  P0,  and  the  density  of  the  fluid,  pj.  From 

this  set  of  eight  independent  variables,  several  dimensionless  parameters  can 
be  formed  that  relate  the  properties  of  the  membrane  to  the  conditions  in  the 
fluid.  The  parameters  include  the  dimensionless  inertia,  M”,  spring  constant, 
K *,  and  membrane  tension,  T*.  The  parameters  are  defined  as  follows: 


M"  = 


IC  = 


T*  = 


m 

pfRo 
KRo 
A  P 
T 

RoAP' 


(19) 

(20) 
(21) 


M*  is  the  ratio  of  the  mass  per  unit  area  of  the  membrane  to  an  equivalent 
mass  per  unit  area  of  the  fluid  based  on  thickness  Ra.  K*  and  Tm  are  the 
ratios  of  the  spring  and  tension  terms,  respectively,  in  the  membrane  model 
to  the  pressure  driving  the  collapse.  In  order  to  obtain  relevant  comparisons 
between  th  e  finite  element  and  finite  difference  results,  it  is  necessary  to 
identify  finite  element  equivalents  for  the  dimensionless  parameters  M*,  Km, 
and  T*. 


The  finite  element  idealization  of  the  spring-backed  membrane  is  shown 
schematically  in  Fig.  5.  There  are  20  elements  in  the  radial  direction  and 
the  membrane  is  one-element  thick.  The  length  of  each  element  is  the  same 
as  that  of  the  corresponding  flow  panels,  i.e.,  0.125  R0.  Symmetry  conditions 
require  that  all  nodes  at  r  =  0  be  constrained  from  moving  in  the  radial  di¬ 
rection.  The  nodes  on  the  edge  of  the  membrane  at  r  =  Rm  are  constrained 
from  moving  vertically.  All  nodes  along  the  base  of  the  backing  elements, 
z  —  —  (<m  +  h),  are  fixed  in  both  the  radial  and  vertical  directions. 


A  very  thin  finite  element  membrane  is  required  to  minimize  bending 
stiffness  and  mass  distribution  effects  that  are  not  present  in  the  finite  differ¬ 
ence  model.  However,  finite  element  aspect  ratio  considerations  require  that 


the  element  thickness  be  no  less  than  one  fourth  of  the  in-plane  dimension. 
Thus,  a  membrane  thickness,  tm ,  of  =  0.125  R0/ 4  =  R0/ 32  is  used.  The 
spring  backing  is  modeled  with  two  rows  of  square  elements  and  has  a  total 
thickness  tb.  The  axial  stiffness  per  unit  area  of  the  backing  elements  is  ob¬ 
tained  from  generalized  Hooke’s  law  (Timoshenko  and  Goodier19),  which  in 
this  case  reduces  to  the  following: 

Ii  =  |  (22) 

where  K  is  the  spring  constant  per  unit  area, 

E  is  Young’s  modulus  of  the  specimen,  and 
t  is  the  thickness  of  the  specimen  in  the  direction  of  the 
applied  load. 


The  combined  axial  stiffness  per  unit  area,  Kt,  of  the  membrane  and 
spring  backing  is  obtained  by  treating  the  two  layers  as  a  set  of  springs  in 


series,  i.e., 


Kt  = 


M-1 

.Em  +  Eb\ 


(23) 


where  Em  and  Eb  are  the  Young’s  moduli  of  the  membrane  and  backing 
elements,  respectively.  The  Young’s  modulus  of  the  backing  elements,  Eb, 
is  chosen  so  that  the  axial  stiffness  per  unit  area  of  the  backing  provides  the 
required  spring  constant,  K.  A  Young’s  modulus  for  the  membrane  is  chosen 
such  that  —j1  >>  Consequently,  the  axial  stiffness  of  the  membrane,  ^ 
has  a  negligible  effect  on  the  combined  axial  stiffness  of  the  membrane  and 
backing.  The  combined  axial  stiffness  per  unit  area  of  the  system  is  now 


(24) 


For  the  finite  element  case,  the  mass  per  unit  area,  m,  is 


TTl  —  pm^m  "I"  pb^b  (25) 

where  pm  is  the  density  of  the  membrane  and  pb  is  the  density  of  the  backing. 
The  density  of  the  backing  elements  is  chosen  so  that  pbtb  «  pmtm]  hence, 
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the  mass  per  unit  area  reduces  to  m  =  pmtm •  The  finite  element  equivalents 
for  M *  and  I\  *  are  expressed  as 


M"  -  pmtm 

PjRo 

(26) 

a?  a, 

II 

a?  a, 

ll 

(27) 

The  finite  element  equivalent  for  T’  is  the  same  as  for  the  finite  difference 
model,  i.e.,  Tm  =  p . 

COMPARISON  OF  RESULTS 

In  the  calculations  described  in  this  section,  the  bubble  is  represented 
by  20  panels  of  equal  arc  length.  The  membrane  radius  is  2.5  R0  and  the 
flow-membrane  interface  is  represented  by  20  panels  of  length  0.125  Ra.  The 
rigid  portion  of  the  flow  boundary  surrounding  the  membrane  is  modeled 
with  40  panels  of  increasing  length,  so  that  the  last  node  is  at  r  =  100  R0. 
This  flow  discretization  is  used  in  both  the  finite  difference  and  finite  ele¬ 
ment  calculations.  The  following  values  are  held  constant  in  both  sets  of 
calculations  to  be  discussed  later;  A P  =  1.0,  pf  —  1.0,  and  Za  —  1.5  R0. 
Furthermore,  the  calculations  are  for  the  collapse  phase  of  the  interaction 
only;  hence  R0  =  Rmax  —  1-0.  In  the  finite  element  calculations,  Poisson’s 
ratio  for  the  spring  backi  ng  i/j  and  that  of  the  membrane  vm  are  taken  as  zero. 

Figure  6  contains  plots  of  the  membrane  velocity  at  r  =  0.0  for  M*  = 
Km  =  2.0.  The  results  of  the  finite  difference  and  finite  element  membrane 
models  are  very  similar.  In  both  cases  the  velocity  reaches  a  maximum  at 
roughly  two  thirds  the  overall  collapse  time  and  decreases  rapidly  at  the  end 
of  the  collapse.  In  Fig.  7,  the  pressure  difference  ( Pm  —  P^)  on  the  membrane 
at  r  =  0.0  is  plotted  for  M*  =  Km  =  2.0.  The  pressure  remains  constant  for 
most  of  the  collapse  and  then  increases  very  rapidly  towards  the  very  end 
of  the  collapse.  For  this  reason,  it  is  most  informative  to  plot  the  pressure 
versus  the  vertical  separation  of  the  north  and  south  poles  of  the  bubble. 
Once  again  the  finite  element  membrane  results  are  in  close  agreement  with 
the  results  of  Duncan  and  Zhang’s  model.  The  pressure  on  the  membrane 
increases  gradually  for  most  of  the  collapse  and  then  rises  rapidly  when  the 
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separation  of  the  poles  of  the  bubble  has  reached  about  20  percent  of  its 
original  value. 


The  excellent  agreement  between  the  results  of  the  two  interaction  models 
provides  a  good  indication  that  the  coupled  finite  element/boundary  element 
model  is  working  properly. 


COMPARISON  WITH  EXPERIMENTAL  DATA 

The  most  extensive  set  of  experimental  data  on  the  behavior  of  vapor  cav¬ 
ities  near  compliant  walls,  is  that  published  by  Shima  et  al.9  What  follows 
is  a  description  of  an  attempt  to  simulate  this  experimental  work  numerically. 

DESCRIPTION  OF  EXPERIMENT 

Figure  8  shows  a  schematic  of  the  experimental  set-up.  The  experiments 
were  conducted  in  a  stainless  steel  bubble  chamber  containing  tap  water  at 
atmospheric  pressure  and  room  temperature.  A  composite  surface  was  placed 
in  contact  with  the  free  surface  of  the  water,  and  a  vapor  bubble  was  then 
generated  by  means  of  electric  spark  discharge.  The  growth  and  subsequent 
collapse  of  the  bubble  was  recorded  photographically.  Composite  wall  prop¬ 
erties  such  as  inertia  and  stiffness  were  varied,  and  their  effect  on  collapse 
behavior  was  observed. 

The  composite  surface  was  40  mm  in  diameter  and  consisted  of  a  Nitrile 
rubber  sheet  of  thickness,  tr,  and  a  foam  rubber  backing  of  thickness,  t  j.  The 
foam  rubber  density,  pj,  was  0.024  g/cm3,  and  the  Nitrile  rubber  density,  pr, 
was  1.39  g/cm3.  A  series  of  uniaxial  compression  tests  was  performed  to 
determine  the  axial  stiffness  of  the  individual  layers.  It  was  observed  that 
for  5  mm  <  f/  <20  mm,  the  spring  constant,  A'/,  for  the  foam  rubber 
was  3.0  kg/cm.  For  the  Nitrile  rubber,  the  value  Kr  =  8.3  kg/cm  was 
reported  for  a  specimen  with  thickness  tr  =  5  mm.  A  set  of  dimensionless 
parameters  was  used  to  describe  the  composite  surface.  The  parameters  were 
the  dimensionless  surface  inertia,  Mm,  and  dimensionless  surface  stiffness,  Km, 
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defined  as  follows: 


(28) 


AT  = 


IC  = 


Prtr  d"  Pftf 
Pi  Rmax 

KtRm  ax 
P  -  p 

*  OO  1  O 


(29) 


where  Kt  is  the  combined  stiffness  of  the  foam  and  Nitrile  rubber  layers 
(Kt  =  ^  ),  Poo  is  the  pressure  in  the  water  at  infinity,  P0  is  the 

pressure  inside  the  bubble  (vapor  pressure  of  the  water),  pi  is  the  den¬ 
sity  of  the  water,  and  Rmax  is  the  maximum  bubble  radius.  In  this  case 
(Poo  —  P<>)  =  1.0159  and  Rma x  —  3.5  mm.  The  thickness  of  the  foam  rubber 
was  held  constant  (tj  =  20  mm).  Variation  of  the  composite  wall  stiffness 
and  inertia  was  obtained  by  varying  the  thickness  of  the  Nitrile  rubber  sheet 
between  0.3  mm  and  5.0  mm. 


NUMERICAL  SIMULATION 


A  schematic  showing  the  numerical  idealization  of  the  experimental  set¬ 
up  is  given  in  Fig.  9.  The  problem  is  axisymmetric  about  the  z-axis.  Thus 
the  problem  can  be  defined  in  the  r-z  coordinate  plane  (0  =  0).  All  lengths 
are  normalized  with  respect  to  the  maximum  bubble  radius,  Rmax -  The  pres¬ 
sure  inside  the  bubble,  P3,  is  assumed  to  equal  the  vapor  pressure  of  the  fluid 
throughout  the  calculation,  i.e.,  Pg  =  Pv  =  const.  The  ambient  fluid  pres¬ 
sure,  Poo,  also  remains  constant  during  the  calculations.  Gravity  effects  have 
not  been  considered;  thus,  the  vapor  bubble  is  shown  above  the  composite 
surface.  The  bubble  is  initially  spherical  with  radius  R0  and  center  at  z  =  Z0. 
Twenty  six  panels  are  used  to  model  the  bubble.  The  interface  between  the 
flow  and  the  surface  of  the  Nitrile  rubber  lies  in  the  r-B  plane  at  z  =  0.  Flow 
along  this  surface  is  modeled  with  40  panels  of  equal  length.  The  flow  bound¬ 
ary  at  z  —  0  is  modeled  as  a  rigid  wall  for  r  >  Rc,  where  Rc  is  the  radius  of 
the  compliant  surface.  The  rigid  portion  of  the  boundary  is  modeled  with  40 
panels  of  increasing  length  so  that  the  last  node  is  at  r  =  100  Rmax ,  z  =  0. 
The  composite  wall  is  modeled  using  finite  elements.  There  are  40  eleme  nts 
in  the  radial  direction,  4  elements  through  the  thickness  of  the  Nitrile  rubber 
layer,  and  8  elements  through  the  foam  rubber  layer.  (A  finite  element  mesh 
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convergence  study  was  performed  for  the  case  of  tT  =  1.0  Rmax  =  1.4). 


It  was  concluded  from  this  effort  that  the  mesh  density  is  sufficient  to 
model  this  thickness  of  Nitrile  rubber.)  Nitrile  and  foam  rubber  densities  are 
as  s  tated  in  the  previous  section.  Young’s  moduli  for  the  two  materials  are 
obtained  from  the  spring  constant  data  for  the  individual  layers.  Applying 
generalized  Hooke’s  law  (Timoshenko  and  Goodier19)  to  the  loading  and 
boundary  conditions  of  the  uniaxial  compression  tests,  yields  the  following 
relation  between  Young’s  modulus  and  axial  stiffness 

a 

t 

where  K  is  the  spring  constant  for  the  Nitrile  or  foam  rubber, 

A  is  the  cross  sectional  area  of  the  rubber  specimen, 
t  is  the  thickness  of  the  rubber,  and 
E  is  the  Young’s  modulus  of  the  material. 


Solving  the  above  relation  for  E  yields  Young’ s  modulus  of  0.330  kPa  and 
0.120  kPa  for  the  Nitrile  and  foam  rubbers,  respectively.  Poisson’s  ratio  for 
the  two  materials  was  not  provided.  Nitrile  rubber  is  nearly  incompress¬ 
ible,  and  a  Poisson’s  ratio,  vT,  of  0.49  is  routinely  used  for  the  treatment 
of  rubber  materials  in  finite  element  analysis.  Conversely,  the  foam  rubber 
is  quite  compressible,  and  the  Poisson’s  ratio  for  the  material  is  probably 
not  a  constant.  A  review  of  data  for  other  foam  rubber  materials  indicates 
a  range  of  Poisson’s  ratio  between  0.24  and  0.36.  Lacking  further  informa¬ 
tion  on  this  particular  material,  a  Poisson’s  ratio  of  vj  —  0.30  is  assumed. 
It  is  important  to  note  that  these  material  properties  were  obtained  from 
static  stiffness  measurements.  However,  strain  rates  cannot  be  expected  to 
remain  near  quasi-static  values  throughout  the  entire  fluid-structure  inter¬ 
action.  Poisson’s  ratio  and  Young’s  modulus  for  many  rubber  and  foam 
materials  are  sensitive  to  the  rate  at  which  the  load  is  applied  (Hunston,  Yu, 
and  Bullman20).  Thus,  the  statically  measured  values  for  Young’s  modulus 
and  Poisson’s  ratio  may  not  accurately  characterize  the  composite  wall  ma¬ 
terials  throughout  the  entire  time  domain  of  interest. 
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COMPARISON  OF  RESULTS 


An  interesting  parameter  describing  the  collapse  of  the  bubble  is  the  col¬ 
lapse  height,  zc.  Collapse  height  is  the  height  above  the  composite  surface 
where  the  north  and  south  poles  of  the  bubble  meet  at  the  end  of  the  col¬ 
lapse.  In  the  published  experimental  data,  the  collapse  height  is  presented 
in  plots  of  zc  versus  the  dimensionless  surface  inertia  M".  Variation  in  M ' 
is  achieved  by  changing  the  thickness  of  the  Nitrile  rubber  layer  while  hold¬ 
ing  the  thickness  of  the  foam  rubber  layer  constant.  Comparative  data  are 
obtained  by  performing  a  series  of  calculations  in  which  the  thickness  of  the 
Nitrile  rubber  elements  in  the  finite  element  model  is  varied.  Note  that  as 
the  thickness  of  the  Nitrile  layer  is  increased,  the  number  of  elements  in  the 
structural  model  increases  also,  while  mesh  density  remains  more  or  less  in¬ 
variant.  The  results  of  the  numerical  study  along  with  the  experimental  data 
for  initial  standoffs  Za  =  1.14  Rmax,  Z0  =  1.43  Rmax  and  Za  =  1.71  Rmax 
are  presented  in  Figs.  10  through  12,  respectively.  For  the  case  of  Za  =  1.14 
Rmax  (Fig.  10),  the  calculated  values  of  collapse  height  are  shifted  towards 
higher  M‘  for  M*  <  1.75.  The  slopes  of  the  calculated  and  experimental 
curves  differ  for  low  M".  However,  at  higher  values  of  M*  the  asymptotic 
behavior  of  the  two  curves  seems  to  be  similar.  For  the  case  of  Za  =  1.43 
Rmax  (Fig.  11),  there  is  good  agreement  between  the  experimental  results 
and  the  present  calculations.  For  Mm  <  1.75,  the  shapes  of  the  two  curves 
are  about  the  same,  with  the  calculations  shifted  to  slightly  larger  M* .  For 
Af*  >  1.75,  the  calculated  collapse  height  drops  off  more  rapidly  than  for  the 
experimental  data.  This  behavior  could  be  an  indication  that  the  finite  ele¬ 
ment  program  is  underpredicting  the  response  of  the  composite  surface  and 
causing  the  interaction  calculation  to  tend  toward  the  rigid  wall  case  more 
rapidly  than  it  should  as  M*  is  increased.  The  collapse  height  predictions 
for  Z0  =  1.71  Rmax  (Fig.  12)  are  slightly  higher  than  the  experimentally 
observed  values;  however,  the  shapes  of  the  two  curves  are  similar. 

There  are  several  possible  sources  of  discrepancy  between  the  experimen¬ 
tal  results  and  the  calculations.  In  the  calculations,  the  pressure  inside  the 
bubble  was  held  constant  at  the  fluid  vapor  pressure.  In  the  experiment, 
there  was  evidence  of  burning  of  electrode  material  during  the  spark  gen¬ 
eration  process.  The  compressibility  of  the  gaseous  products  resulting  from 

18 


electrode  burning  is  likely  to  have  affected  the  experimental  results.  The 
inadequacy  of  the  linear  elastic  material  model  to  represent  the  Nitrile  rub¬ 
ber  is  likely  to  be  a  major  source  of  error  in  the  calculations.  A  strain-rate 
sensitive  nonlinear  elastic  representation  of  the  rubber  properties  is  a  neces¬ 
sity  in  these  calculations.  Finally,  the  linearized  boundary  condition  on  the 
fluid-structure  interface  is  likely  to  have  influenced  the  predictions  for  low 
Af*.  In  these  cases,  the  Nitrile  surface  attains  a  significant  normal  velocity, 
but  this  was  not  considered  in  the  linearized  Bernoulli  equation. 

INTERACTION  WITH  A  COMPLIANT  SPHERE 
MODEL  DESCRIPTION 

The  interaction  calculations  discussed  thus  far  have  involved  planar  fluid- 
structure  interfaces,  constant  bubble  internal  pressure,  and  no  gravitational 
effects.  Attention  is  now  focused  on  a  more  complex  problem  involving  a 
doubly-  curved  fluid-structure  interface,  namely,  the  interaction  between  a 
growing  and  collapsing  bubble  and  a  submerged  spherical  structure.  The 
problem  is  shown  schematically  in  Fig.  13.  The  problem  is  rotationally  sym¬ 
metric  and  is  defined  in  cylindrical  coordinates  (r-z)  at  0  =  0.  The  z-axis 
pierces  the  north  and  south  poles  of  both  bubble  and  sphere.  The  spherical 
structure  has  outer  radius  R,,  wall  thickness  t,  and  is  centered  at  r  =  0, 
z  —  —Rt.  By  requiring  the  spherical  structure  to  be  thin  walled,  >  50, 
thin  shell  theory  can  be  U3ed  to  validate  the  structural  portion  of  the  model. 
As  in  the  previous  calculations,  the  bubble  is  initially  spherical  with  radius 
Ro  and  center  at  r  =  0,  z  =  Z0  (the  center  of  the  bubble  is  initially  at 
depth  d  below  the  fluid  free  surface)  and  the  pressure  inside  the  bubble  is 
initially  Pg  <>.  Bubble  internal  pressure  is  made  to  obey  the  poly  tropic  law, 
Pa  —  Pgo  (^)  ,  where  Vo  is  the  initial  bubble  volume  and  V  is  the  instan¬ 
taneous  bubble  volume.  The  bubble  is  modeled  with  30  flow  panels  and  the 
”wet”  shell  interface  is  modeled  with  60  panels.  The  structural  finite  element 
model  has  60  meridional  elements  and  one  element  through  the  shell  thick¬ 
ness  (Fig.  14).  The  structural  nodes  on  the  axis  of  symmetry  are  prevented 
from  moving  radially  but  are  free  to  move  in  the  vertical  direction.  Thus  the 
spherical  structure  can  undergo  rigid  body  vertical  motion,  as  well  as  local 
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deformation. 


There  are  10  independent  variables  in  this  interaction  problem.  For  the 
flow,  the  variables  are  the  init.’  d  pressure  differential,  A P  =  Pgo  —  Poo:  fluid 
density,  pj;  initial  position  of  the  bubble,  Z0\  initial  radius  of  the  bubble, 
R0\  and  maximum  bubble  radius,  Rmax ■  ^or  the  spherical  structure,  the 
variables  are  the  outer  radius,  R,;  wall  thickness,  t3;  Young’s  modulus,  Es; 
Poisson’s  ratio,  v3 ;  and  density  of  the  structural  material,  ps.  An  interesting 
parameter  in  this  problem  is  the  collapse  height,  zc,  which  is  the  height  above 
the  north  pole  of  the  sphere  where  the  north  and  south  poles  of  the  bubble 
meet  at  the  end  of  the  collapse.  Collapse  height  is  a  function  of  the  10 
independent  variables  previously  discussed,  but  can  be  expressed  in  terms  of 
seven  dimensionless  parameters 

is  =  f  (Jk-  _^L_  is  _iA_  p.  (&'-*$)  K  B  \ 

Zo  *  \  Rmax  ’  Rmax  ’  R*  ’  -Rmax  ’  PfRmax*  ’  Rm.x  AP ’  R3max  AP  )  1  j 


where  Ri  is  the  inner  radius  of  the  shell  (Ri  =  Ra  —  ta), 

K  is  its  extensional  rigidity  (Kraus21), A'  =  ,  and 

B  is  its  bending  rigidity  (Timoshenko  and  Woinowsky-Krieger22) 


B  =  - 
12(1 


For  ease  of  reference,  the  last  three  shell  terms  on  the  right-hand  side  of 
Eq.  31,  are  defined  as  M*,  A'*,  and  B‘,  respectively.  M *  is  the  mass  ratio 
of  the  sphere  to  the  displaced  fluid  at  Rmax •  The  dimensionless  membrane 
stiffness,  Km,  is  the  ratio  of  the  extensional  stiffness  to  the  pressure  driving 
the  collapse.  The  dimensionless  flexural  stiffness,  Bm,  relates  bending  stiffness 
to  the  pressure  driving  the  collapse.  In  this  study,  the  outer  radius  of  the 
spherical  structure  is  a  multiple  of  Rmax Ra  =  nfCoi •  Thus,  B *  can  be 
expressed  in  terms  of  the  radius  of  the  shell 

„._n3  E  f ta\ 3  1 

12(1  —  v2)  \RaJ  A P’ 

It  should  be  noted  that  since  Ra/ta  >  50  and  1  >>  (f,/i?,)3,  the  bending 
stiffness  of  the  sphere  is  relatively  small  and  its  influence  here  is  negligible. 


RESULTS 


The  time  scale  for  the  flow  in  this  problem  is  Ta  =  Rmax\Jj In  the  cal¬ 
culations  that  are  described  in  this  section,  all  times  and  lengths  are  nondi- 
mensionalized  by  T0  and  Rmax ,  respectively.  Furthermore,  pj,  P^,  and  Rmax 
are  taken  as  unity  in  their  respective  units;  y  =  1.25;  Z0  =  1.5  Rmax  and 
the  calculations  begin  with  the  bubble  radius  Ra  =  0.4  Rmai.  The  initial 
bubble  internal  pressure  is  Pg0  =  1.5  P^.  The  bubble’s  initial  conditions 
have  been  chosen  to  be  consistent  with  the  detonation  of  one  half  pound  of 
TNT  at  a  depth  of  300  ft  (see  Cole3  for  a  detailed  discussion  of  the  gas  globe 
that  results  from  the  detonation  of  TNT).  The  structure  has  radius  R,  =  4 
Rmax  and  wall  thickness,  t,  --  0.075  Rmax  so  that  R3/ta  =  53.33.  For  now, 
gravitational  and  depth  effects  are  not  included  in  the  calculations.  The  am¬ 
bient  hydrostatic  pressure,  P*,,  is  applied  statically  to  the  structure  at  the 
start  of  the  transient  calculation.  The  resulting  hydrostatic  displacements 
are  part  of  the  initial  conditions  for  the  structure  in  the  subsequent  dynamic 
calculation.  The  deformation  of  the  shell  due  to  this  static  preload  must  be 
small  since  the  boundary  conditions  are  satisfied  at  the  undisturbed  position 
of  the  fluid-  structure  interface.  This  condition  dictates  a  fairly  stiff  shell 
(high  A'");  hence,  Young’s  modulus  for  the  sphere  is  fixed  at  E ,  =  10,000.0 
Poo-  The  resulting  value  of  Km  is  then  4286,  and  the  static  deflection  due  to 
Poo  is  6  =  f,/10.0  for  Poisson’s  ratio  vt  =  0.30. 

We  now  consider  the  effect  of  Mm  on  the  collapse  of  the  bubble  while 
holding  fixed  all  other  dimensionless  parameters  and  again  neglecting  gravi¬ 
tational  effects.  Variation  in  M*  is  achieved  by  changing  the  density  of  the 
structural  material.  Figure  15  contains  plots  of  the  height  of  the  north  and 
south  poles  of  the  bubble  as  a  function  of  time.  The  plots  are  for  four  dif¬ 
ferent  values  of  M'  and  for  the  case  in  which  the  sphere  is  rigid  and  fixed  in 
space.  The  height  at  which  the  poles  of  the  bubble  meet  is  indicated  by  a 
short  dashed  line.  In  the  rigid  sphere  case,  the  collapse  height  is  zc  =  0.966 
Rmax-,  and  the  total  time  for  the  growth  and  collapse  is  tc  =  1.856  T0.  For  a 
flat  rigid  wall  with  Ra  =  0.4  Rmax  and  Za  =  1.5  Rmax,  the  collapse  height  is 
zc  =  0.85  Rmax,  and  the  collapse  time  is  tc  =  2.024  Ta.  Thus,  adding  curva¬ 
ture  to  the  fluid-structure  interface  significantly  alters  the  dynamic  behavior 
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of  the  bubble.  When  the  sphere  is  compliant  but  very  massive  (A/*  =  226), 
the  collapse  time  is  tc  —  2.104  Ta,  and  the  collapse  height  is  zc  =  0.865  Rmax. 
Reducing  the  dimensionless  mass  of  the  sphere  to  M *  =  92  yields  tc  and  zc 
values  of  2.104  T0,  and  1.16  Rmax ,  respectively.  As  M *  is  decreased  further, 
the  collapse  height  increases  monotonically,  and  the  collapse  time  decreases. 
At  these  lower  values  of  M’,  the  poles  of  the  bubble  do  not  meet  at  the 
end  of  the  collapse  calculation.  This  is  an  indication  that  the  collapse  has 
become  almost  spherical.  (The  collapse  height  for  a  bubble  collapsing  spher¬ 
ically  in  an  infinite  fluid  is  zc  =  Z0.)  It  should  be  noted  that  in  these  cases  of 
near  spherical  collapse,  the  calculation  becomes  unstable  near  the  end  of  the 
collapse;  hence,  the  final  position  of  the  bubble  is  a  little  more  speculative 
than  in  the  other  calculations.  As  the  mass  of  the  sphere  is  reduced  below 
M *  =  92,  the  separation  of  the  poles  at  the  end  of  the  collapse  begins  to  de¬ 
crease,  and  the  projected  collapse  height  is  greater  than  the  initial  standoff, 
indicating  redirection  of  the  bubble  motion. 

Profiles  of  the  bubble  at  various  times  during  the  collapse  near  the  rigid 
sphere  and  several  compliant  spheres  of  various  M*  values  are  given  in  Fig. 
16.  In  the  profiles  for  the  rigid  sphere  case,  a  reentrant  jet  directed  toward 
the  sphere  can  be  seen  near  the  end  of  the  collapse.  For  M*  =  226,  the  later 
stages  of  collapse  are  also  characterized  by  the  formation  of  a  reentrant  jet 
directed  towards  the  sphere.  As  mentioned  above,  the  collapse  height  is  ac¬ 
tually  smaller  than  the  collapse  height  for  the  rigid  sphere  case.  When  M*  is 
reduced  to  76  (Fig.  16),  the  collapse  has  been  modified  such  that  all  collapse 
profiles  are  nearly  spherical.  The  collapse  height  in  this  case,  zc  =  1.34  Rmax , 
is  beginning  to  approach  the  value  for  spherical  collapse  in  an  infinite  fluid, 
zCoo  =  1.5  Rmax •  When  M*  is  reduced  to  56  (Fig.  16),  the  collapse  height, 
zc  =  1.53  Rmax,  is  higher  than  the  initial  standoff,  Z0  =  1.5  Rmax-  In  this 
case,  deviation  from  sphericity  is  apparent.  The  shape  of  the  bubble  near 
the  south  pole  suggests  that  a  redirected  jet  is  beginning  to  form  late  in  the 
collapse. 

The  effect  of  M*  on  the  response  of  the  spherical  structure  is  illustrated 
in  Fig.  17,  which  contains  plots  of  dimensionless  velocity  of  the  north  pole, 
in  the  shell  normal  direction  (positive  outward),  versus  time.  As  expected, 
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shell  peak  velocity  increases  as  the  sphere  mass  is  reduced.  Figure  18  con¬ 
tains  similar  plots  for  the  south  pole.  In  all  cases,  the  south  pole  velocity 
is  lower  than  that  of  the  north  pole,  indicating  that  interaction  between  the 
bubble  and  sphere  is  driven  by  local  deformation  of  the  surface  of  the  sphere 
near  the  bubble.  North  pole  displacement  of  the  shell  is  shown  in  Fig.  19 
for  several  values  of  M*.  Peak  displacements  increase  as  structure  mass  is 
decreased.  For  M‘  =  56,  the  peak  structural  displacement  (0.066  Rma x) 
is  nearly  90  percent  of  the  shell  thickness,  t3.  Local  displacements  of  the 
sphere  in  execess  of  one  shell  thickness  violate  the  underlying  assumptions  of 
small  displacements  of  the  fluid-structure  interface,  and  hence,  the  numeri¬ 
cal  model  fails.  For  this  reason,  calculations  for  M*  <  56  were  not  performed. 

The  effect  of  gravity  is  now  included  in  the  calculations  by  expressing  the 
gravity  term  in  Bernoulli’s  equation  (Eq.  2)  in  terms  of  the  dimensionless 
Froude  number,  i.e.,  Eq.  (2)  becomes 

W  =  iW)’  +  ^+<i^  (33) 

where  all  quantities  are  dimensionless  and  the  Froude  number  is  defined  by 
Fr  =  Pool PjgRmax  —  d/Rmax  where  d  is  the  initial  depth  of  the  centroid  of 
the  bubble.  Fig.  20  contains  plots  of  the  height  of  the  north  and  south  poles 
of  the  bubble  vs  time  for  M*  =  226  and  various  values  of  Fr.  When  FT  =  206, 
i.e.,  the  bubble  is  initially  206  Rmax  below  the  fluid  surface,  the  effect  of  grav¬ 
ity  is  negligible,  and  the  collapse  time  and  collapse  height  are  essentially  the 
same  as  for  the  case  when  gravity  was  not  considered  (d  =  206  Rmax  is  con¬ 
sistent  with  the  intitial  conditions  used  in  the  previous  set  of  calculations). 
As  the  initial  bubble  depth  is  reduced  to  100  Rmax ,  gravitational  effects  have 
increased  the  collapse  height  to  0.95  Rmax •  For  the  case  of  d  =  10  Rmaxi  grav¬ 
itational  effects  have  become  more  important  than  interaction  forces,  and  the 
collapse  height  becomes  greater  than  2.0  Rmax-  Figure  21  contains  collapse 
profiles  for  M*  =  226  at  various  depths.  For  the  cases  of  d  =  206  Rmax  and 
d  =  100  Rmaxi  significant  target  attraction  and  target-  directed  jetting  occur. 
However,  when  d  is  reduced  to  10  Rmaxi  the  collapse  is  driven  by  gravita¬ 
tional  effects,  and  the  bubble  migration  and  jetting  are  away  from  the  target. 
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SUMMARY 


For  a  given  initial  standoff,  collapse  behavior  near  a  rigid  spherical  target 
is  significantly  different  from  that  near  a  flat  rigid  wall.  When  the  sphere  is 
very  compliant,  the  collapse  problem  is  sensitive  to  the  mass  of  the  spherical 
structure.  When  the  target  mass  is  high,  a  reentrant  jet  forms  late  in  the 
collapse  and  is  directed  toward  the  structure.  As  the  shell  mass  is  decreased, 
the  collapse  becomes  spherical  with  no  jet  formation.  A  jet,  directed  away 
from  the  target,  begins  to  form  when  the  structure  mass  is  low.  At  very 
low  structural  mass,  significant  local  shell  deformation  occurs,  and  the  small 
displacement  assumptions  for  the  fluid-structure  interface  are  violated.  When 
gravitational  forces  are  considered,  the  collapse  is  driven  by  interaction  effects 
at  depths  greater  than  100  Rmax-  The  collapse  is  driven  by  gravitational 
effects  for  shallow  depths. 

CONCLUSION 

In  this  study  a  boundary  element  flow  algorithm  was  coupled  with  a  finite 
element  structural  analyzer  to  perform  interactive  calculations  of  the  growth 
and  collapse  of  an  explosion  bubble  near  a  compliant  structure.  The  valid¬ 
ity  of  the  program  was  established  through  direct  comparisons  with  a  code 
previously  developed  by  Duncan  and  Zhang,10,11  in  which  the  structure  was 
a  spring-backed  membrane.  There  was  good  agreement  between  the  results 
of  the  two  programs. 

Calculations  of  the  collapse  of  a  vapor  bubble  next  to  a  composite  wall 
were  performed.  The  calculated  collapse  heights  were  in  good  agreement 
with  published  experimental  data. 

A  parametric  study  of  the  interaction  between  an  explosion  bubble  and  a 
thin  walled  spherical  structure  was  also  performed.  The  mass  of  the  spherical 
structure  was  found  to  be  a  critical  parameter  in  the  interaction.  When  the 
mass  of  the  structure  was  high,  the  collapse  of  the  bubble  was  characterized 
by  the  formation  of  a  reentrant  jet  directed  toward  the  structure.  As  the 
mass  of  the  spherical  structure  was  decreased,  the  collapse  became  spherical 
and  no  jet  formed.  A  jet,  directed  away  from  the  spherical  structure,  began 


to  form  when  the  mass  of  the  structure  was  low.  Significant  shell  deforma¬ 
tion,  enough  to  violate  small  displacement  assumptions,  was  associated  with 
very  low  structure  mass.  The  ratio  of  the  depth  of  submergence  to  bubble 
maximum  radius  was  also  found  to  be  a  critical  parameter  in  the  collapse 
problem.  When  this  ratio  is  large  (greater  than  100),  the  collapse  is  driven 
by  interaction  forces.  However,  for  shallow  submergence,  gravitational  effects 
become  more  important  than  interaction  forces.  For  a  given  initial  standoff, 
the  collapse  of  a  bubble  near  a  rigid  sphere  was  found  to  be  significantly 
different  from  the  case  of  collapse  next  to  a  rigid  wall. 

A  promising  technique  has  been  synthesized  for  numerical  prediction  of 
coupled  explosion  bubble-structure  interaction.  Extension  to  a  full  three- 
dimensional  model  will  make  possible  submarine  and  surface  ship  whipping 
calculations  with  completely  interactive  loading,  target  attraction,  and  buoy¬ 
ancy  effects.  This  method,  in  combination  with  a  shock  wave  fluid-structure 
decoupling  scheme,  will  eventually  make  possible  the  prediction  of  structural 
response  to  both  shock  wave  and  hydrodynamic  loading  from  a  pulsating 
bubble  experiencing  gravity  migration  as  well  as  target /free- surface  attrac¬ 
tion  or  repulsion. 
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FLUID  FREE-SURFACE 


Fig.  i.  Schematic  showing  the  coordinate  system  and  the  initial  position 
of  the  bubble  and  the  compliant  spherical  shell. 
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Fig.  3.  Nodes  and  panels  on  the  surface  of  the  bubble  and  on  the 
fluid-structure  interface. 
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INITIAL  POSITION  OF  CAVITY 


Flfl.4.  Schematic  showing  spring  backed  membrane  and  initial  configuration 
of  the  vapor  bubble. 
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Fig.  5.  Finite  element  idealization  of  spring  backed  membrane. 
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Fig.  6.  Dimensionless  normal  velocity  of  membrane  at  r  *  0.0  for  M*  =  K*  =  2.0. 
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Fig.  7.  Dimensionless  pressure  Pffl/AP  on  membrane  at  r  =  0.0  versus  separation 
of  north  and  south  poles  of  the  bubble  (Znp-z^/Rmax  f°r  M*  =  K*  =  2.0. 
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40  mm 


BRASS  BACKING  PLATE 


Fig.  8.  Schematic  showing  experimental  set-up  for  compliant  coating 
experiments  of  Shima,  Tomita,  Gibson,  and  Blake. 


z 


ALL  NODES  FIXED  ALONG  BASE 
OF  THE  FORM  RUBBER 


Fig.  9.  Finite  element  idealization  of  compliant  coating. 
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Fig.  10.  Comparison  of  results  of  present  calculations  and  data  from  experiments  of 

Shima,  Tomita,  Gibson  and  Blake.  Plot  of  Zc/Rmax  versus  M*  for  Zo  =  1 .14  R^. 


Fig.  11.  Comparison  of  results  of  present  calculations  and  data  from  experiments  of 

Shima,  Tomita,  Gibson  and  Blake.  Plot  of  Zc/Rmax  versus  M*  for  Zo  =  1 .43  R^. 
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Fig.  13.  Schematic  showing  coordinate  system,  spherical  shell,  and  initial 
position  of  the  explosion  bubble. 
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OMENS  IONLESS  HEIGHT 


DIMENSIONLESS  TIME 


Fig.  15.  Height  of  north  and  south  poles  of  the  bubble  versus  time. 
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3.00 


GROWTH  AND  COLLAPSE  PROFILES  FOR  RIGID  SHELL 


RADIAL  COORDINATE 


Fig.  16.  Bubble  collapse  profiles  at  various  times  for  collapse  near  rigid  shell  and 
compliant  shells  of  M*  =  226,  76,  and  56. 
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Fig.  16.  (Continued) 
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Fig.  17.  Velocity  (positive  outward  into  fluid)  of  the  north  pole  of  the  shell  versus  time. 
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Fig.  18.  Velocity  (positive  outward  into  fluid)  of  the  south  pole  of  the  shell  versus  time. 
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DIMENSIONLESS  TIME 


Effect  of  gravity  on  the  motion  of  the  north  and  south  poles  of  the  bubble. 
Plot  of  height  of  north  and  south  poles  of  the  bubble  versus  time  for  shells 
of  M*  *  226,  at  depths,  d  =  206  Rmax,  100  R^,  10  R^. 
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RADIAL  COORDINATE 

Fig.  21.  Bubble  collapse  profiles  at  various  times  near  shell  of  M*  =  226,  at  depths, 

d  *  206  Rfnaxi  100  Rmax,  10  Rmax- 
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